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1  Project  Summary 

Although  ignored  by  most  molecular  biologists  who  commonly  characterize  gene  ex¬ 
pression  levels,  the  reality  is  that  the  level  of  expression  of  the  same  gene  can  vary 
enormously  from  one  cell  to  another  within  a  genetically-identical  cell  population.  It 
has  been  shown  that  these  fluctuations  are  not  simply  a  by-product  of  the  regulatory 
process  and  that  they  can  contribute  significantly  to  the  control  of  cell  function.  For 
example,  some  organisms  utilize  fluctuations  to  introduce  diversity  into  a  population, 
as  occurs  during  the  lysis-lysogeny  transition  in  A  phage,  while  in  others  stability 
against  fluctuations  is  essential  in  gene  regulatory  cascades  controlling  processes  such 
as  differentiation.  This  project  developed  improved  theoretical  and  simulation  tech¬ 
niques  that  lead  to  reliable  predictions  of  complex  cell  behaviors.  Moreover,  this 
project  developed  quantitative  models  of  gene  regulation  that  correctly  incorporate 
stochastic  fluctuations  that  are  inherent  in  all  gene  networks  and  currently  limit  the 
usefulness  of  existing  modeling  approaches.  Computational  models  were  tested  and 
refined  through  experimental  studies  using  engineered  gene  circuits.  Through  this 
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integrated  approach,  a  hierarchical  simulation  approach  was  developed  and  incorpo¬ 
rated  into  BioSPICE  in  order  to  enhance  the  user’s  ability  to  describe,  predict  and 
control  complex  gene  networks  in  living  cells. 


2  Significant  Accomplishments 

The  significant  accomplishments  for  the  project  are  summarized  below: 

•  GRASS  v2.0.  Gene  Regulatory  Anaylsis  and  Stochastic  Simulation,  was  devel¬ 
oped,  successfully  integrated  and  released  as  part  of  BioSPICE  v2.0. 

•  SDEsolver,  a  general  purpose  biochemical  reaction  simulator,  was  developed, 
successfully  integrated  and  released  as  part  of  BioSPICE  v2.0. 

•  BioNetS,  a  biochemical  network  simulator,  was  developed,  successfully  inte¬ 
grated  and  released  as  part  of  BioSPICE  v3.0.  BioNetS  handles  discrete,  con¬ 
tinuous,  and  mixed  models.  It  also  provides  a  semi-implicit  method  for  stiff 
systems,  and  is  optimized  for  speed  and  validated  for  accuracy. 

•  Models  for  sources  of  noise  (transcription  and  translation)  in  gene  expression 
systems  (single-gene  systems  and  cascade  networks  of  genes)  were  developed 
and  experimentally  validated  using  engineered  gene  networks. 

•  Models  of  synthetic  engineered  promoters  were  developed  and  experimentally 
validated  using  engineered  gene  circuits. 

•  A  synthetic  engineered  promoter,  which  could  be  independently  activated  and 
repressed  by  two  different  proteins,  was  constructed  and  used  to  validate  the 
developed  models. 

•  Gene  expression  systems  (single-gene  systems  and  cascade  networks  of  genes) 
were  constructed  and  used  to  validate  models  for  sources  of  noise  (transscription 
and  translation)  in  gene  expression. 

•  A  “noise  generator”  that  can  be  used  to  produce  gene  regulatory  signals  with 
variable  signal-to-noise  ratios  was  developed. 

•  Prototype  programmable  cells  in  which  engineered  regulatory  modules  were 
integrated  with  the  cell’s  natural  circuitry  were  developed. 
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3  Software  for  Stochastic  Modeling  of  Biochemical 
Networks 

3.1  Background 

Mathematical  modeling  of  complex  biological  networks  has  a  lengthy  history.  In  the 
past,  the  standard  approach  for  modeling  these  systems  has  been  to  derive  ordinary 
differential  equations  (ODEs)  based  on  the  law  of  mass  action  for  the  concentra¬ 
tions  of  the  biochemical  species  involved  in  the  network.  Experimental  studies  have 
demonstrated,  however,  that  stochastic  effects  can  be  significant  in  cellular  reactions, 
particularly  in  the  case  of  transcriptional  regulation,  where  generally  there  are  two 
copies  of  each  gene  and  the  number  of  messenger  RNA  (mRNA)  molecules  can  be 
small.  A  number  of  recent  experimental  and  modeling  studies  have  addressed  the  role 
of  fluctuations  in  gene  expression.  Many  modeling  studies  have  employed  the  well- 
established  Gillespie  Monte  Carlo  algorithm  or  one  of  its  more  recent  variants.  These 
algorithms  offer  an  exact  solution  to  the  stochastic  evolution  of  chemical  systems, 
but  they  are  computationally  very  expensive.  A  much  more  efficient  approach  is  to 
approximate  the  species  as  continuous  variables  and  formulate  the  problem  in  terms 
of  stochastic  differential  equations  (SDEs),  often  referred  to  as  chemical  Langevin 
equations.  This  approximation  works  remarkably  well  for  many  cases,  even  when  the 
number  of  particles  involved  is  as  small  as  10,  and  the  resulting  simulations  can  run 
orders  of  magnitude  more  quickly  than  the  discrete  Monte  Carlo  approach.  In  other 
cases,  when  some  or  all  of  the  particle  numbers  are  very  small,  the  system  may  need  to 
be  modeled  using  the  discrete  approach,  or  a  hybrid  method  in  which  some  species  are 
treated  discretely  while  others  are  evolved  using  the  continuum  approximation.  With 
the  increasing  interest  in  formulating  accurate  models  of  large  biochemical  networks, 
there  is  a  need  for  reliable  software  packages  that  correctly  incorporate  stochastic 
effects,  yet  are  fast  enough  to  simulate  large  interconnected  sets  of  reacting  species 
(as  found,  for  example,  in  signaling  cascades  or  genetic  regulatory  networks).  The 
Biochemical  NETwork  Stochastic  Simulator,  “BioNetS,”  was  developed  to  meet  this 
need.  BioNetS  is  capable  of  performing  full  discrete  simulations  using  an  efficient 
implementation  of  the  Gillespie  algorithm.  It  is  also  able  to  set  up  and  solve  the 
chemical  Langevin  equations,  which  are  a  good  approximation  to  the  discrete  dy¬ 
namics  in  the  limit  of  large  abundances.  Finally,  BioNetS  can  handle  hybrid  models 
in  which  chemical  species  that  are  present  in  low  abundances  are  treated  discretely, 
whereas  those  present  at  high  abundances  are  handled  continuously.  Thus,  the  user 
can  pick  the  simulation  method  that  is  best  suited  to  their  needs.  All  aspects  of  the 
software  are  highly  optimized  for  efficiency. 

In  the  Implementation  section,  the  mathematical  background  for  the  Gillespie 
method,  chemical  Langevin  equations  and  hybrid  models  is  presented,  along  with  a 
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discussion  of  the  numerical  algorithms  used  in  BioNetS.  Under  Results  and  Discussion, 
a  brief  introduction  to  BioNetS  is  provided  along  with  several  examples.  The  examples 
serve  two  purposes:  (1)  to  illustrate  how  to  use  the  software,  and  (2)  to  verify  its 
efficiency  and  accuracy. 

3.2  Implementation 

The  mathematical  methodology  on  which  BioNetS  is  built  is  developed  first. 

3.2.1  Discrete  Reactions  and  the  Gillespie  Algorithm 

BioNetS  makes  use  of  elementary  reactions  (zeroth,  first  and  second  order).  The 
following  examples  illustrates  each  type  of  reaction: 


0 

SL 

Y 

A 

(1) 

A 

B 

(2) 

Y 

A  +  B 

AJ3 

(3) 

k4 

V 

kQ 

V  +  V 

(4) 

In  the  above  reactions,  the  calligraphic  letters  denote  a  single  molecule  of  a  chemical 
species.  The  number  of  molecules  of  a  particular  species  in  the  system  at  time  t 
is  denoted  with  uppercase  letters  (e.g.,  A(t),  R(t),  and  All  the  rate 

constants,  7,  S,  and  ki-ke,  have  units  of  per  time.  Eq.  1  represents  a  process  in  which 
a  molecule  A  is  produced  when  the  reaction  proceeds  in  the  forward  direction  and 
is  degraded  in  the  reverse  direction.  In  the  forward  direction,  the  reaction  is  zeroth 
order  and  proceeds  with  an  average  rate  of  7.  In  the  backward  direction,  the  reaction 
is  first  order,  and  the  average  rate  of  degradation  is  5A{t).  The  forward  reaction  in 
Eq.  2  represents  a  process  in  which  chemical  species  A  is  converted  to  species  B.  In 
this  case,  A  and  B  might  represent  two  different  conformations  of  the  same  molecule. 
In  Eq.  2  both  the  forward  and  backward  reactions  are  first  order  because  the  reaction 
rates  are  proportional  to  the  respective  concentrations.  The  forward  reaction  given 
in  Eq.  3  is  a  second-order  reaction  in  which  an  A  molecule  and  a  B  molecule  come 
together  to  form  the  complex  AJ3.  The  average  rate  for  the  reaction  is  kiA{t)B{t). 
The  backward  reaction  is  a  first-order  reaction  in  which  A-B  dissociates  at  an  average 
rate  of  k2A_B{t).  In  Eq.  4  the  forward  reaction  produces  a  molecule  V.  The  difference 
between  this  reaction  and  the  forward  reaction  in  Eq.  1  is  that  the  average  rate  is 
kA^{t).  This  leads  to  exponential  growth  of  V{t).  This  reaction  is  particularly  useful 
if  V (t)  is  interpreted  as  the  cell  volume.  In  the  backward  reaction,  two  V  molecules 
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come  together  and  degrade  one  of  the  V  molecules.  The  average  rate  for  this  reaction 
is  —  1).  The  V{t)  —  1  term  arises  because  two  of  V{t)  molecules  must  be 

chosen  to  react.  This  type  of  term  also  arises  in  reactions  that  produce  homodimers. 
This  reaction  eventually  stops  the  exponential  growth  of  V.  The  net  effect  of  these 
two  reactions  is  to  produce  logistic  growth.  The  total  average  reaction  rate  for  the 
set  of  reactions  given  in  Eqs.  1-4  is 


A^(^)  —  7  T  SA(t)  +  kiA(t)  +  k2B{t)  + 

ksA{t)B{t)  +  kiA_B{t)  +  k^V{t)  +  keV{t){V{t)  -  1) 

4 

=  ^(F,  +  Bt)  (5) 

i=l 

where  Fj  and  Bi  are  the  average  forward  and  backward  rates,  respectively,  for  the  ith 
reaction. 

For  the  rest  of  this  section,  it  is  assumed  that  the  volume  of  the  cell  is  not  changing 
and  only  Eqs.  1-3  are  considered.  In  the  examples,  a  case  in  which  the  volume  is 
changing  is  considered.  If  A{t)^  B{t)  and  A_B{t)  are  present  in  large  numbers,  then  the 
law  of  mass  action  can  be  applied  to  derive  equations  that  govern  the  concentrations 
[A]  =  A{t)/V,  [B]  =  B{t)jV  and  [AB]  =  A_B{t)/V,  where  V  is  the  cell  volume. 
These  equations  are 


dt 

4B] 

dt 

d[AB] 

dt 


—  {k'^[B]  ki  <5)[^]  -|-  A;4[^F]  -|-  A;2[F]  -|-  ^y' 
-{k'^[A]  +  k2)[B]  +  ki[A]  +  ki[AB] 
k',[A][B]  -  k,[AB] 


(6) 

(7) 

(8) 


The  primed  rate  constants  indicate  that  they  have  been  appropriately  scaled  by  the 
volume  (i.e,  k'^  =  k^V  and  7'  =  7/E),  and,  therefore,  have  units  of  either  per  time 
per  concentration  or  concentration  per  time.  Note  that  to  convert  to  units  of  molar, 
one  also  has  to  appropriately  scale  the  rate  constants  by  Avagadro’s  number.  Eqs.  6  - 
8  represent  a  macroscopic  description  of  the  process,  because  they  ignore  fluctuations 
in  the  concentration  that  arise  from  the  stochastic  nature  of  chemical  reactions. 

In  general,  A(t),  B{t)  and  A_B{t)  are  random  variables  that  take  on  any  nonneg¬ 
ative  integer  value.  The  Gillespie  algorithm  can  be  used  to  generate  sample  paths  of 
the  process.  This  algorithm  assumes  that  the  random  time  ATj,  between  the  ith  and 
i  -I-  1  reaction,  is  exponentially  distributed.  For  the  simple  example  given  by  Eqs.  1 
-  3,  the  mean  waiting  time  between  reactions,  which  characterizes  the  exponential 
distribution,  is  =  7  +  SA{ti)  -|-  kiA{ti)  -|-  k2B{ti)  -|-  kzA{ti)B{ti)  -|-  kiA_B{ti), 
where  ti  is  the  time  at  which  the  ith  reaction  occurred.  Therefore,  ti+i  —  ti  +  ^Ti. 
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Once  the  time  at  which  the  next  reaction  occurrs  has  been  determined,  the  following 
probabilities  are  used  to  determine  which  reaction  occurred: 


Pr[0^A]  = 

7 

PATi 

(9) 

Pr[A^0]  = 

SA{ti) 

PATi 

(10) 

Pt[A^B]  = 

kiA{ti) 

PATi 

(11) 

Pr[B^A\  = 

k2B{ti) 

PATi 

(12) 

Pt[A  +  B^  A^]  = 

kiA{ti)B{ti) 

PATi 

(13) 

Pt[AJ3^  A  +  B]  = 

kiA_B{ti) 

PATi 

(14) 

Once  the  reaction  has  been  determined,  the  chemical  species  are  updated  accordingly. 

As  discussed  in  the  Numerical  Methods  section,  BioNetS  uses  an  efficient  implemen¬ 
tation  of  the  Gillespie  algorithm. 

Another  description  of  discrete  stochastic  processes  is  achieved  through  use  of  the 
master  equation  that  governs  how  the  probabilities  of  the  various  random  variables 
in  the  process  evolve  in  time.  Let  Pa,b,ab{t)  =  Pr[A(t)  =  a,B{t)  =  b,A_B{t)  =  a_6], 
then  Pa,b,aj){t)  satisfies  the  master  equation 

—  [7  +  (5  +  ki)a  +  k2b  +  k^ab  +  k4aJ)]pa,b,a.b  +  'yPa-l,b,aJ>  +  +  l)Pa-|-l,6,a.6 

-\-ki{a  +  l)Pa+l,6-l,a.6  +  +  l)Pa-l,6+l,a.6 

+k3{a  +  1)(6  +  l)Pa+l,6+l,a.6-l  +  ki{aJ)  +  l)Pa-l,6-l,a.6-|-l  (15) 

The  master  equation  is  the  starting  point  for  deriving  various  approximate  schemes  for 
describing  the  system.  In  the  next  section,  an  approximate  scheme  that  is  valid  in  the 
limit  of  large,  but  finite  molecule  numbers,  is  discussed.  The  simplest  approximation 
scheme  is  achieved  by  considering  the  first  moments  of  the  process.  Over  bars  will  be 
used  to  denote  averaging.  For  example,  A{t)  =  J2abajb'^Pa,b,a-b{'t)-  Eq-  15  can  be  used 
to  derive  equations  that  govern  the  time  evolution  of  all  the  first  moments.  Because 
of  the  second-order  reaction  in  Eq.  3,  the  equations  for  the  means  are  coupled  to 
the  second  moments.  In  fact,  the  nth  moment  equations  contain  terms  that  involve 
the  n  -I-  1  moments.  Thus,  there  is  no  closure  to  the  system.  The  simplest  closure 
scheme  is  to  assume  that  all  moments  factorize  (e.g.,  A^  =  A  ).  This  represents 
the  macroscopic  limit  in  which  fluctuations  are  ignored.  In  this  limit,  Eqs.  6-8  are 
recovered  from  the  master  equation. 


dpa,b,aJ) 

dt 
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The  Diffusion  Limit  and  the  Chemical  Langevin  Equations 

The  general  form  of  the  master  equation  for  a  system  consisting  of  N  chemical  species 
and  M  reactions  is 

,  M  MM 

+  Bi)pn  +  ^  FiPn-Si  +  BiPn+Si  (16) 

i=l  i=l  i=l 

where  n  is  a  N-dimensional  vector  of  species  numbers,  F*  and  are  the  backward 
and  forward  rates  for  the  ith  reaction,  and  the  vectors  5i  contain  the  stoichiometric 
constants  for  the  ith  reaction.  For  the  simple  model  given  by  Eqs.  1  -  3,  iV  =  3, 
M  =  3,  and  Pnit)  =  Pr[^(t)  =  =  n2,  and^_F(t)  =  ua].  The  forward  and 

backward  rates  are  Fi  =  7,  =  Srii,  F2  =  kiUi,  B2  =  k2n2.  Fa  =  kzn\n2-,  and 

Bz  =  kJ^nz.  and  the  5,  vectors  are  the  rows  of  the  stoichiometric  matrix 

1  0  0  \ 

-110 
-1-11/ 


The  (i,  j)  element  in  the  above  matrix  represents  the  change  in  the  jth  chemical 
species  when  the  ith  reaction  proceeds  in  the  forward  direction. 

If  the  molecule  numbers  are  large  as  compared  to  1,  then  the  master  equation 
Eq.  16  can  be  approximated  by  the  continuous  process 


5p(n,  t) 


dt 


where 


M 


‘^j,i {Fj  Bi) 


M 


^i,j^i,k{Fi  Bi) 


(18) 


(19) 

(20) 


This  result  can  be  derived  in  several  ways.  One  method  is  to  note  that  Eq.  15 
represents  a  second-order  finite  differencing  of  Eq.  18,  with  a  grid  size  of  1.  Another 
method  is  to  make  use  of  the  shift  operator 

o  00  -jj  o 

f{n  +  k)  =  exp{k—)f{n)  =  ^  ^^/(^)  (21) 

where  f{n)  is  an  arbitrary  smooth  function  and  for  our  purposes  k  is  an  integer.  If 
the  shift  operator  is  used  in  Eq.  15,  the  diffusion  limit  is  achieved  when  the  Taylor 
series  expansion  given  in  Eq.  21  is  truncated  at  j  =  2. 
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Sample  paths  consistent  with  Eq.  18  can  be  generated  using  the  following  set  of 
SDEs 

M 

=  >1,(N)  +  Y,  Ai„VFUNy+a(N)«.t(t)  (22) 

fe=l 

where  the  Wk{t)  are  independent  Gaussian  white  noise  processes.  These  equations 
are  often  referred  to  as  the  chemical  Langevin  equations.  For  Eqs.  1-3,  the  explicit 
form  of  the  SDEs  are 

=  —{k\  +  k^B  +  6)A.  +  k2B  +  k^A^B  -1-  q  +  s/ 5A  +  'ywi(t) 

-\/kiA  +  k2Bw2{t)  -  ^/ksAB  +  kiA_Bwz{t)  (23) 

=  —{k2  +  kzA)B  +  k\A  +  kiA_B  +  sJkiA  +  k2Bw2{t) 

—  \/kzAB  +  kiA^wz  (t)  (24) 

=  — k/^A^B  T  k^AB  T  ■sJk^AB  T  kiA_Bw2{t^  (25) 

BioNetS  generates  numerical  solutions  to  the  SDEs  given  by  Eq.  22  using  either  an 
explicit  or  semi-implicit  Euler  method.  The  form  of  these  methods  is 

Nj{t  +  At)  =  Nj{t)  +  At^^-(N(t  +  eAt)) 

M 

+ VAt  ^  Ak,j  y/Fk{N{t))  +  Bk{N{t))Zk{t)  (26) 

k=l 

where  e  =  0  for  the  explicit  method  and  e  =  1  for  the  semi-implicit  method  and  the 
Zk{t)  are  independent  standard  normal  random  variables.  The  advantage  of  using  the 
chemical  Langevin  equations  is  that  in  the  appropriate  parameter  regime,  numerical 
solutions  to  the  set  of  SDEs  given  by  Eq.  22  can  be  generated  much  more  efficiently 
than  using  the  Gillespie  algorithm.  This  point  is  expanded  upon  in  the  Results 
section.  Higher-order  numerical  algorithms  for  SDEs  are  available,  but  the  noise 
structure  of  the  chemical  Langevin  equations  makes  these  schemes  very  cumbersome 
to  implement.  In  the  Results  section,  the  Euler  methods  given  by  Eq.  26  are  verified 
to  be  sufficient  to  produce  reliable  results.  Note  that  the  A  matrix  is  generally  sparse, 
and  BioNetS  takes  advantage  of  this  sparseness  to  optimize  the  efficiency  of  the  two 
Euler  methods  (see  Numerical  Methods,  below). 

3.2.2  Hybrid  Schemes 

It  is  often  desirable  to  allow  some  of  the  chemical  species  to  be  treated  as  continuous 
random  variables  and  some  to  be  treated  discretely.  This  is  particularly  true  for  the 
case  of  transcriptional  regulation  by  transcription  factors.  In  this  situation  there  can 
be  as  few  as  one  DNA/transcription  factor  binding  site  and  mRNA  abundances  can  be 


dA 

dt 

dB 

dt 

dA_B 

dt 
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as  small  as  10  or  fewer.  In  contrast,  protein  abundances  can  be  in  the  thousands.  The 
technical  difficulty  with  implementing  hybrid  schemes  that  include  both  discrete  and 
continuous  random  variables  is  that  the  Gillespie  method  requires  constant  transition 
rates  between  reactions.  This  may  not  be  the  case,  if  some  of  the  chemical  species  are 
evolving  continuously  in  time.  BioNetS  overcomes  this  problem  in  one  of  two  ways. 

Let  Nd  <  N  he  the  number  of  discrete  chemical  species  and  Md  <  M  the  number 
of  reactions  that  produce  a  change  in  one  of  the  Nd  chemical  species.  The  overall 
reaction  rate  at  time  tj  for  the  discrete  set  of  chemical  species  is 

Md 

#■«,  =  (27) 

i=l 

If  the  time  step  At  for  the  SDEs  is  small  enough  such  that 


Pt  =  HtjNt  <  e  «  1  (28) 

then  Pt  is  approximately  the  probability  of  a  transition  in  At.  In  the  above  equation 
e  is  a  user-specified  tolerance.  The  probability  of  two  discrete  transitions  in  At 
is  proportional  to  (At)^.  Choosing  e  <  0.1,  which  means  the  probability  of  two 
reactions  in  At  is  less  than  0.01,  generally  produces  good  results.  However,  this 
should  be  verified  on  a  case-by-case  basis.  At  each  time  step,  BioNetS  checks  to 
verify  that  Ineq.  28  is  satisfied  for  the  specified  e.  If  so,  a  uniform  random  number 
R  is  generated  and  compared  against  pt-  li  R  <  pt,  then  a  transition  occurred  and 
the  conditional  probability  R/pt  is  used  to  determine  which  of  the  discrete  transitions 
occurred.  If  pt  >  e,  then  the  discrete  reactions  determine  the  fastest  time  scale  in  the 
system.  In  this  case  the  Gillespie  algorithm  is  used  to  update  the  discrete  reactions, 
and  the  random  time  step  Ntj  is  used  to  update  the  SDEs. 

3.2.3  Numerical  Methods 

BioNetS  generates  code  that  is  tailored  to  efficiently  simulate  biochemical  reactions. 
The  optimization  techniques  used  by  BioNetS  allows  the  software  to  simulate  large 
systems  in  reasonable  times  without  requiring  high-end  computational  hardware. 
Techniques  used  to  optimize  the  Gillespie  method  are: 

•  For  the  discrete  variables,  the  program  uses  data  structures  that  allow  only  the 
chemical  species  and  reaction  rates  that  are  affected  by  the  current  reaction  to 
be  updated. 

•  A  bisection  search  is  used  to  determine  which  reaction  occurred. 

The  code  has  both  an  explicit  and  a  semi-implicit  solver,  for  simulating  the  chem¬ 
ical  Langevin  equations.  The  user  specifies  at  runtime  which  method  to  use.  By 
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default  the  semi-implicit  solver  will  be  used.  The  semi-implicit  solver  uses  Newton’s 
method  to  solve  the  implicit  equations,  and  for  that  the  program  needs  to  compute 
the  Jacobian  and  solve  a  linear  system  at  each  iteration.  For  updating  the  chemical 
Langevin  equations  and  hybrid  models  optimization  techniques  include: 

•  The  sparse  nature  of  the  stoichiometric  matrix  is  used  to  efficiently  store  and 
perform  matrix  operations. 

•  After  every  reaction,  only  the  species  and  reaction  rates  affected  by  that  reaction 
are  updated. 

•  The  Jacobian  is  sparse,  and  the  code  takes  full  advantage  of  this  fact.  The 
program  solves  and  factorizes  the  Jacobian  using  sparse  methods.  Before  the 
code  generation,  BioNetS  computes  the  entries  in  the  Jacobian  symbolically 
and  finds  a  permutation  that  decreases  the  number  of  fill-ins  during  the  LU 
factorization.  As  a  result,  no  zero  entries  are  saved,  and  the  sparse  structure  is 
fully  exploited.  The  sparse  structure  is  then  used  in  the  LU  solve.  In  the  code, 
no  pivots  are  visible,  and  no  if-statements  are  left. 

3.3  Results  and  Discussion 

In  this  section,  several  examples  which  serve  as  illustrations  of  how  to  use  BioNetS 
and  test  the  accuracy  and  efficiency  of  the  numerical  methods  are  presented.  One 
particular  concern  is  the  accuracy  of  the  Euler  methods.  While  these  methods  are 
only  of  order  a/^,  it  is  shown  that  when  the  approximations  that  lead  to  the  chemical 
Langevin  equations  are  valid,  the  difference  between  the  numerical  solutions  of  the 
SDEs  and  the  exact  discrete  Gillespie  method  are  negligible.  Currently,  the  graphical 
user  interface  to  BioNetS  runs  on  the  Macintosh  OS  X  operating  system,  though  the 
software  will  generate  portable  C/C++  code  that  can  be  compiled  and  run  in  any 
computing  environment.  The  following  examples  illustrate  the  way  in  which  models 
are  entered  and  run  in  BioNetS.  More  detailed  documentation  is  available  with  the 
software  package. 

3.3.1  Dimerization 

To  begin,  consider  a  simple  system  that  consists  of  the  following  two  reactions: 


0 

A 

M 

(29) 

<5m 

M  +  M 

V 

(30) 

V 

0 

(31) 
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In  this  system,  monomer  molecules  A4  are  produced  at  an  average  rate  7  and  degraded 
at  an  average  rate  Two  monomers  can  then  bind  to  form  a  dimer  molecule 

V.  The  average  forward  and  backward  rates  for  this  reaction  are  —  1) 

and  k2D{t),  respectively.  The  dimers  are  degraded  at  a  rate  5d-  Two  cases  will  be 
treated.  In  the  first  case  the  cell  volume  is  assumed  to  be  constant,  and  in  the  second 
case  the  cell  is  allowed  to  grow  and  divide.  To  model  cell  growth,  the  cell  volume  I4 
is  treated  as  a  random  variable  I4  =  OiV ,  where  1/  is  a  non-negative  discrete  random 
variable  and  a  represents  a  unit  of  volume.  The  random  variable  V  is  governed  by 
the  reaction 

v4v  +  V  (32) 

The  above  reaction  causes  V  to  grow  exponentially  fast  with  an  average  rate  of 
ks-  Note  that  logistic  growth  is  produced  when  the  backward  reaction  in  Eq.  32  is 
included. 

3.3.2  Constant  Volume 

Start  by  considering  the  simple  case  in  which  the  volume  of  the  cell  remains  constant. 
To  use  BioNetS  follow  these  steps.  Copy  BioNetS  onto  your  machine,  and  double  click 
to  launch.  Help  is  included  as  part  of  the  program,  and  accessed  from  the  Help  menu. 
The  Help  document  will  walk  you  through  all  the  steps  needed  to  enter  reactions  and 
run  the  simulator. 

The  user  interface  asks  you  to  enter  the  reaction  and  corresponding  rate  constants 
in  the  top  part  of  the  script  window.  In  the  bottom  part  of  the  script  window,  you  can 
toggle  between  panels.  The  Species  panel  allows  the  user  to  specify  how  the  simulator 
treats  each  chemical  species,  discrete  or  continuous.  The  Constants  panel  lists  the 
order  in  which  the  rate  constants  are  referenced.  The  Output  panel  allows  the  user 
to  specify  the  ouput  type.  There  are  two  ways  to  generate  program  output,  either 
binary  or  ASCII.  Binary  output  is  based  on  MATLAB  binary  files,  so  it  is  possible  to 
drive  the  program  with  MATLAB  and  use  MATLAB ’s  plotting  routines  to  view  the 
output.  It  is  also  possible  to  generate  time  series  and  histograms  of  the  data  from 
within  BioNetS.  Using  ASCII  files  for  I/O  allows  the  simulator  to  be  run  through  shell 
scripts.  The  Executable  panel  allows  the  user  to  generate  either  an  executable  file 
or  source  code.  BioNetS  generates  portable  C/C++  code  that  can  be  compiled  and 
run  in  any  computing  environment.  BioNetS  can  directly  compile  the  C/C++  code. 
However,  this  requires  the  Developer  tools,  included  on  all  recent  Apple  machines 
and  available  directly  from  http://developer.apple.com  for  free.  The  compiled  code 
can  then  be  run  from  within  BioNetS.  The  Comments  panel  is  available  for  the  user 
to  enter  descriptive  comments  about  the  model. 

To  run  BioNetS  as  a  BioSpice  agent,  you  need  to  move  the  source  directory  onto 
a  OAA-supported  system.  Once  there,  open  up  the  MakeOAA  file  and  specify  the 
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locations  of  your  oaalib  folder.  Then  just  type  “make  -f  MakeOAA”  (without  the 
quotes)  to  create  the  agent. 

Simluations  indicated  that  the  agreement  between  the  two  different  methods  is 
very  good.  These  findings  indicate  that  the  chemical  Langevin  equations  are  accu¬ 
rately  capturing  the  dynamics  and  steady-state  behavior  of  the  discrete  system. 

3.3.3  Cell  Growth  and  Division 

In  this  section,  it  is  described  how  cell  growth  and  division  can  be  modeled  using 
BioNetS.  It  is  assumed  that  the  cell  is  experiencing  exponential  growth  up  until  the 
time  it  divides.  As  discussed  above,  the  cell  volume  I4  is  treated  as  a  random  variable. 
In  this  model  cell  division  occurs  when  Vc  exceeds  a  threshold  value  When  cell 
division  occurs  the  volume  is  halved,  and  the  proteins  are  randomly  divided  between 
the  two  cells  using  a  binomial  distribution.  Only  one  of  the  daughter  cells  is  tracked. 
Because  second-order  reactions  require  two  molecules  to  collide,  the  rate  constants  for 
these  reactions  should  scale  like  ki  =  k[jVc.  It  is  also  assumed  that  the  production 
rate  of  monomers  scales  as  7  =  q'W.  This  is  a  reasonable  assumption,  because  as  the 
cell  grows  the  transcription  and  translation  machinery  increases.  These  assumptions 
produce  the  following  rate  equations  for  the  concentrations 

-2k'^[Mf  +  2k2[D]  +  7'  -  {5^  +  k^)[M]  (33) 

k'^[Mf-k2[D]-{k^  +  5d)[D]  (34) 

kzV,  (35) 

The  terms  in  Eqs.  33  and  34  that  involve  kz  arise  because  of  dilution  due  to  cell 
growth.  The  same  parameter  values  as  in  the  constant  volume  case  are  used  except 
5m  =  I  and  5d  =  0.  The  cell  growth  rate  kz  =  0.02  and  the  scale  factor  for  the 
volume,  a,  is  equal  to  1.  With  these  choices  of  parameter  values,  we  expect  the 
average  behavior  of  this  system  to  be  similar  to  that  of  the  constant  volume  case. 
Simulations  for  this  simple  example  indicated  that  the  main  effect  of  volume  growth 
is  to  act  as  an  additional  noise  source  and  increase  the  variability  of  the  distributions. 

3.3.4  A  Chemical  Oscillator 

BioNetS  is  next  used  to  simulate  a  two-gene  system  that  has  been  studied  in  the 
literature.  In  this  system,  the  protein  A  coded  for  by  gene  a  acts  as  an  activator  for 
gene  a  and  gene  r,  by  binding  to  the  promoter  regions,  Va  and  Vr:  of  the  respective 
gene.  This  increases  the  rate  of  mRN Aa  and  mRN Ar  production  by  a  factor  and 
Or,  respectively.  The  protein  K  acts  as  a  repressor  for  both  genes  by  binding  to  A  to 
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form  the  inactive  complex  AJZ.  All  gene  products,  mRNA  and  protein,  are  actively 
degraded.  However,  the  heterodimer  AJZ  protects  the  V.  subunit  from  degradation. 
The  system  consists  of  9  chemical  species  and  the  following  14  biochemical  reactions: 


Va 

N 

Va  +  mRNAa 

(36) 

Vu-A 

agfcl 

VaNl  +  mRNAa 

(37) 

Vr 

Vr  +  mRNAr 

(38) 

Vr^ 

ark2 

Vr-A  +  mRNAr 

(39) 

mTZMAa 

mRNAa  +  A 

(40) 

mRNAr 

mRNAr  +  R 

(41) 

A  -\-  R- 

A-R 

(42) 

kQ 

Va+  A 

Va-A 

(43) 

ks 

Rr  T  A 

k^ 

Vr-A 

(44) 

kio 

A 

fci} 

0 

(45) 

R 

fcig 

0 

(46) 

mRNAa 

fci? 

0 

(47) 

mRNAr 

ki^ 

0 

(48) 

AJl 

ki^ 

R 

(49) 

An  interesting  feature  of  the  system  is  that  it  is  capable  of  producing  sustained 
oscillations. 

The  chemical  species  Va,  Vr,  Vr-A,  and  Vr-A  are  binary  random  variables:  they 
can  only  take  on  the  values  0  or  1.  Therefore,  these  species  cannot  be  approximated 
as  continuous  random  variables.  All  the  other  chemical  species  appear  in  sufficient 
quantities  to  justify  the  continuum  approximation.  The  hybrid  model  was  run  using 
the  semi-implicit  Euler  method,  and  for  these  parameter  values,  runs  three  times 
faster  than  full  model.  Visually,  the  agreement  between  the  two  methods  appears 
good.  To  test  the  accuracy  of  the  Euler  method,  BioNetS  was  used  to  construct  2-D 
histograms  of  R  versus  mRNAr. 

Simulations  showed  excellent  agreement  between  the  discrete  and  hybrid  models. 
This  indicates  that  the  hybrid  model  is  accurately  sampling  the  steady-state  distribu¬ 
tion.  To  verify  that  the  hybrid  model  faithfully  captures  the  dynamics  of  the  system, 
the  power  spectra  of  both  models  were  computed.  Again,  excellent  agreement  was 
found  between  the  discrete  and  hybrid  models. 
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3.3.5  An  Engineered  Promoter  System 

Genetic  regulatory  networks  consist  of  sets  of  genes  whose  levels  of  expression  in  a 
cell  are  interdependent.  This  dependence  arises  through  the  action  of  transcription 
factors,  proteins  which  bind  to  operator  sites  on  the  DNA  strand  and  influence  the 
rate  at  which  a  gene  product  is  generated.  Once  bound,  these  regulatory  proteins 
affect  the  binding  affinity  of  RNA  polymerase,  an  enzyme  that  binds  to  promoter  se¬ 
quences  in  the  DNA  and  initiates  transcription  of  messenger  RNA  (mRNA)  strands, 
which  are  subsequently  translated  into  proteins.  These  proteins  may  themselves  act 
as  transcription  factors,  influencing  their  own  rates  of  expression  or  those  of  other 
gene  products  and  thus  forming  networks  of  connected  genes.  Using  standard  tech¬ 
niques  in  modern  molecular  biology,  it  is  possible  to  design  novel  systems  of  promoter- 
gene  pairs,  such  that  virtually  any  desired  regulatory  network  architecture  may  be 
instantiated;  such  networks  are  often  called  “synthetic  gene  networks.”  Recent  im¬ 
plementations  have  included  direct  negative  and  positive  feedback,  a  bistable  switch, 
an  oscillator,  an  intercellular  communication  system,  and  a  bimodal  self-activating 
system. 

In  this  example,  BioNetS  is  used  to  implement  a  model  of  a  simple,  open-loop 
network  based  around  a  novel  engineered  promoter,  which  was  designed  and  con¬ 
structed  by  the  project  team.  The  promoter,  called  OjjO/ac,  combines  the  Oiac:  OrI^ 
and  Or2  operator  sites,  so  that  it  is  repressed  by  the  lac  repressor  protein  (Lad) 
and  activated  by  the  lambda  repressor  protein  (Cl).  Experiments  were  conducted  in 
which  the  promoter,  along  with  other  sites  to  produce  the  activator  and  repressor 
proteins,  is  integrated  into  a  high  copy  number  plasmid  and  inserted  into  a  strain 
of  Escherichia  coli.  The  promoter’s  activity  is  observed  using  a  fluorescent  reporter. 
Green  Fluorescent  Protein  (GFP).  The  goal  here  is  to  provide  a  reasonably  complex 
test  case  to  evaluate  the  performance  of  BioNetS. 

The  processes  to  be  captured  by  the  model  are:  transcription  and  degradation 
of  mRNA  strands;  translation  of  mRNA  into  protein;  degradation  of  protein;  forma¬ 
tion  of  protein  multimers  (dimers  in  the  case  of  Cl,  tetramers  in  the  case  of  Lad); 
Lad  binding  to  isopropyl- /3-D-thiogalactopyranoside  (IPTG),  a  chemical  inducer  that 
reduces  Lad’s  binding  affinity  for  Oiac]  and  protein-DNA  binding  at  the  OnOiac  pro¬ 
moter’s  operator  sites.  We  define  the  following  chemical  species:  G,  GFP;  Mg,  mRNA 
coding  for  GFP;  X,  Cl  monomer;  X2,  Cl  dimer;  mRNA  coding  for  Cl;  D^,  the 
arabinose-inducible  pBAD  promoter  site  producing  Cl;  Y,  Lad  monomer;  1^2,  Lad 
dimer;  Yi,  Lad  tetramer;  /q,  IPTG  (present  in  massive  excess  and  thus  taken  to  be 
constant);  Yj,  Lad  tetramer  bound  to  IPTG;  My,  mRNA  coding  for  Lad;  and  Dy,  the 
PitetOl  site  constitutively  producing  Lad.  In  addition  to  these,  species  Dq  through 
Ds  are  defined,  representing  the  various  permutations  of  proteins  bound  to  the  three 
operator  sites  in  the  OjjO/ac  promoter.  There  are  twelve  combinatorial  possibilities. 
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but  three  of  them  are  eliminated  on  the  basis  that  Cl  {X2)  binding  Or2  but  not  OrI 
is  unlikely,  because  of  the  low  binding  affinitity  of  Cl  for  Or2  compared  to  Oijl.  This 
reflects  the  regulatory  effect  of  the  proteins;  for  example,  Cl  bound  to  Or2  leads  to 
a  10-fold  increase  in  transcription  rate,  while  Lad  bound  to  Oiac  halts  transcription 
completely  (note  that  we  assume  in  the  event  of  simultaneous  binding  of  activator 
and  repressor,  repression  “wins”  and  transcription  is  halted). 

The  following  irreversible  reactions  represent  the  processes  of  transcription,  trans¬ 
lation,  and  degradation: 


Vo 

Vq  +  My 

(50) 

v^ 

Vl  +  J\4y 

(51) 

V2 

V2  +  My 

(52) 

Vy 

Vy  -\-  My 

(53) 

V, 

V3;  +  Mx 

(54) 

My 

My  +  g 

(55) 

My 

My+y 

(56) 

Xtx 

M-x  +  X 

(57) 

My 

'ymr^a 

0 

(58) 

My 

'ymr^a 

0 

(59) 

•M.X 

'ymr^a 

0 

(60) 

g 

yprot 

0 

(61) 

y 

yprot 

0 

(62) 

X 

yprot 

0 

(63) 

As  in  previous  reactions,  the  caligraphic  letters  represent  individual  molecules  of  each 
species.  All  times  and  rates  are  scaled  by  the  cell  division  time. 

Experimental  measurements  generally  provide  equilibrium  rather  than  rate  con¬ 
stants,  and  thus  when  writing  reversible  reactions  we  use  the  following  notational 
convention:  a  reaction  with  equilibrium  constant  K  has  forward  rate  constant  KR 
and  backward  rate  constant  i?,  where  i?  is  a  scaling  factor  which  sets  the  speed  at 
which  the  reaction  approaches  equilibrium  (three  values  of  i?  -  1,  10,  and  100  - 
are  considered).  Using  this  notation,  protein-protein  binding  is  represented  with  the 
following  set  of  reactions 
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+  (64) 

Ky2 

:^^2  +  :t^2  ^  (65) 

y4  +  io  -  (66) 

A:  +  A:  P  A'a  (67) 

Finally,  protein-DNA  binding  is  given  by: 

2^0  +  ^2  ^  2^1  (68) 

2^1  +  A2  2?2  (69) 

2^o  +  :1^4  ^  2^3  (70) 

2^i  +  :1^4  ^  2^4  (71) 

2^3 +  A2  ^  Vi  (72) 

2^2  +  :1^4  ^  V,  (73) 

2^4 +  A2  P  2^5  (74) 

Vo  +  yi  P  2^6  (75) 

2^1  +  3^/  ^  2^7  (76) 

2^6 +  A2  P  Vj  (77) 

2^2  +  :1^7  ^  2^8  (78) 

r>7  +  ^2  ^  2^8  (79) 


In  all,  the  system  consists  of  21  species,  participating  in  34  reactions.  The  reactions 
are  entered  into  BioNetS  using  the  same  method  described  in  the  previous  examples. 
BioNetS’  ability  to  represent  individual  species  as  either  discrete  or  continuous  is  used 
to  formulate  three  versions  of  the  model:  fully  discrete,  fully  continuous,  and  a  hybrid 
version  in  which  the  DNA  species  Dq  through  are  discrete  while  all  other  species 
are  continuous.  The  value  of  2?,  the  scaling  factor  for  reversible  reactions  is  varied, 
and  all  other  parameters  are  kept  fixed  at  the  following  nondimensionalized  values: 
/3g  =  0.1,  I3y  =  1,  /3^  =  0.5,  /3t  =  10,  Jrnma  =  3.5,  Jprot  =  0.7,  Ky  =  0.01,  Ky2  =  0.1, 
Kyi  =  2x  10-®,  =  0.05,  Ki  =  0.3,  K2  =  2Ki,  K3  =  0.008,  =  1.4  x 

Jo  =  1  X  10®. 

To  evaluate  the  steady-state  probability  distributions  produced  by  the  reaction 
system,  simulations  250000  cell  cycles  in  length  were  used  to  accumulate  histograms 
(a  built-in  feature  of  BioNetS)  of  the  number  of  molecules  of  GFP  (species  G),  for 
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each  of  the  three  versions  of  the  model.  The  resulting  distributions  were  essentially 
identical,  indicating  that  the  continuum  approximations  used  in  the  fully  continuous 
and  hybrid  forms  of  the  model  were  valid.  Not  all  species  in  the  system  are  well 
approximated  as  continuous  variables  The  fully  continuous  model  fluctuates  into  neg¬ 
ative  values,  indicating  that  the  continuum  approximation  has  broken  down.  This 
does  not  significantly  affect  the  distribution  for  GFP  because  the  other,  more  com¬ 
mon  DNA  states  dominate  the  system’s  behavior;  note,  however,  that  if  genomic  DNA 
were  considered  rather  than  a  high  copy  number  plasmid,  one  would  not  be  able  to 
employ  a  fully  continuous  model.  The  hybrid  model,  by  treating  the  DNA  species  as 
continuous,  eliminates  the  fluctuations  into  negative  values.  In  general,  the  appro¬ 
priate  approximations  will  depend  on  both  the  system  and  the  variables  of  interest: 
in  the  present  example,  if  one  were  interested  in  the  behavior  of  the  operator  sites 
themselves,  one  would  not  be  able  to  use  the  fully  continuous  version  of  the  model, 
but  as  a  model  solely  of  GFP  expression  the  approximation  suffices.  Comparisons 
between  types  of  models  should  be  made  to  test  the  underlying  assumptions,  and 
BioNetS  facilitates  this  process. 

Simulations  200  cell  cycles  in  length  were  used  to  test  the  speed  at  which  the  three 
model  versions  ran.  In  each  case,  200  simulations  were  run  using  a  consistent  set  of 
200  different  random  seeds;  all  runs  were  started  with  identical  initial  conditions.  For 
the  fully  continuous  and  hybrid  systems,  the  semi-implicit  scheme  was  numerically 
stable  and  yielded  consistent  histograms  for  all  time  step  sizes  between  dt  =  0.001 
and  dt  =  0.5,  but  the  latter  corresponds  to  just  two  time  points  per  cell  division 
cycle  (recall  that  all  times  are  scaled  by  the  cell  division  time),  and  it  was  chosen 
instead  to  sample  20  points  per  cycle  and  set  dt  =  0.05.  Simulations  showed  that 
the  fully  continuous  method  was  always  fastest,  with  the  degree  of  improvement  over 
the  exact,  fully  discrete  method  depending  strongly  on  the  value  of  i?,  the  scaling 
factor  for  the  reversible  reaction  rates.  For  i?  =  1,  the  fully  continuous  method 
was  only  1.4-fold  faster  than  the  fully  discrete  method,  but  as  R  is  increased  this 
speed  advantage  increases  to  over  4-fold  at  R  =  10,  then  to  over  30-fold  at  R=  100. 
(Note  that  the  speed  advantage  of  the  fully  continuous  over  the  fully  discrete  method 
increases  with  the  abundances  of  the  chemical  species.  Shifting  parameters  to  generate 
higher  protein  numbers  can  yield  cases  in  which  the  continuum  approximation  is 
hundreds  of  times  faster  than  the  discrete  approach;  runs  not  shown  here.)  Use 
of  a  hybrid  discrete/continuous  method  did  not,  for  this  particular  model  system, 
offer  any  speed  gain  over  the  fully  discrete  approach;  the  increased  time  involved  in 
computing  the  Jacobian  for  the  semi-implicit  method  is  more  time-consuming  than 
simply  simulating  the  reactions  directly.  Optimizing  efficiency  requires  testing  various 
potential  approaches,  and  BioNetS  makes  this  a  simple  process. 
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3.4  Conclusions 

BioNetS  was  developed  to  be  a  reliable  tool  for  studying  the  stochastic  dynamics  of 
large  chemical  networks.  The  software  allows  the  user  to  specify  which  of  the  chemi¬ 
cal  species  in  the  network  should  be  treated  as  discrete  random  variables  and  which 
can  be  approximated  as  continuous  random  variables.  The  software  is  highly  opti¬ 
mized  for  speed  and  should  be  be  able  to  simulate  networks  consisting  of  hundreds  of 
chemical  species.  The  accuracy  of  the  numerical  methods  was  verified  by  considering 
several  test  systems  (a  dimerization  reaction,  a  chemical  oscillator,  and  an  engineered 
promoter),  each  of  which  shows  excellent  agreement  between  the  fully  discrete  ver¬ 
sion  and  the  fully  or  partially  continuous  versions.  BioNetS,  by  providing  a  simple, 
user-friendly  interface,  will  allow  biological  experimentalists  to  formulate  biochemical 
reaction  models  of  their  systems  quickly  and  easily,  ideally  increasing  the  number  of 
systems  in  which  direct  comparisons  are  available  between  models  and  experimental 
results.  Clearly,  not  every  possible  biological  system  can  be  captured  in  the  current 
version  of  BioNetS,  and  its  capabilities  will  continue  to  grow  in  the  future. 

3.5  Requirements 

•  Project  name:  BlOchemical  NETwork  Stochastic  Simulator  (BioNetS) 

•  Operating  system: 

-  User  interface:  Macintosh  OS  X,  version  10.2  or  above. 

-  Generated  source  code:  Ability  to  compile  portable  C-|— I-  code.  Makefiles 
included  for  OS  X  and  Linux. 

•  Programming  language:  C-|— 1-. 

•  Other  requirements:  None. 
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